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Abstract 

This paper presents some numerical simulations of the full one-dimensional Maxwell-Lorentz equations 
that describe light propagation in fiber Bragg gratings in order to confirm that the standard nonlinear coupled 
mode equations fail to predict the weakly nonlinear dynamics of the system when dispersive instabilities 
come into play, and that, in this case, the correct slow envelope description of the system requires to consider 
higher order dispersion effects. 
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I. INTRODUCTION 



The nonlinear coupled mode equations (NLCME) are the envelope equations currently 
used to study the weakly nonlinear dynamics of light propagation in fiber Bragg gratings. 
These equations do not include dispersion effects. In this paper we integrate numerically the 
full ID Maxwell-Lorentz equations in a fiber grating in order to show that the dispersion 
effects can be essential in the dynamics of the system and that the correct weakly nonlin- 
ear description of the system has necessarily to include higher order dispersion terms. The 
resulting envelope equations are asymptotically nonuniform in the sense that they include 
terms with different asymptotic order, and this is a standard situation for general extended, 
propagative (i.e., with order one group velocity) pattern forming systems. 

The weakly nonlinear dynamics of resonant light propagation in a Fiber Bragg grating (FBG), 
i.e., optical fiber with a periodic variation of the refractive index along its length, is usually de- 
scribed using the so-called nonlinear coupled mode equations (NLCME) 

At= A+ + iKA- +iA+{a\A+\^ + \A-f), (1) 
At = -A- + iKA+ + iA-{a\A-\^ + \A+\^), (2) 

which prescribe the evolution of the complex envelopes A"^ of the two slowly modulated resonant 
wavetrains that approximately constitute the actual field inside the FBG 

A+ (x, t) e^"+^"* + A~ (x, t)e^^"+'"* + c.c. + ■ ■ • , (3) 

see e.g. The NLCME above, where space, time and the amplitudes have been 

rescaled to reduce the number of parameters, retain the combined effect of the group velocity, 
the coupling induced by the grating and the weakly nonlinear interaction of the wavetrains. This 
formulation, apart from FBG, has been also used to describe the evolution of quasi-onedimensional 



Bose-Einstein condensates in optical lattices |^ |70 and, in general, the NLCME are commonly 
regarded as the normal form for the weakly nonlinear dynamics of any extended, propagative 
system without dissipation and with a weak spatial periodic structure. 

In a recent paper [8] one of the authors showed that the NLCME ©-© fail to predict the 
dynamics of the system when dispersive instabilities (that cannot be detected using the NLCME 
formulation) come into play, and that, for both signs of the dispersion coefficient, there are always 
stable solutions according to the NLCME that are dispersively unstable. In order to correctly 
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describe the weakly nonlinear evolution of the system, the effect of higher order dispersion has to 
be retained and the appropriate amplitude equations are the following dispersive nonlinear coupled 
mode equations (NLCMEd) 

At = At + inA- + iA+{a\A+\^ + \A-\^) + i£A+ , (4) 
A^ = -A- + ■mA+ + iA-{a\A-\^ + \A+\^) + isA-^. (5) 



This dispersive system was introduced and analyzed in detail in [|8|], but we think it is convenient 
to briefly remind here some of the results obtained in that paper: 

1. The scaling of the NLCMEd is the same of the NLCME: the characteristic length scale is 
the slow scale that results from the balance of the advection term with the small effect of 
the grating, the characteristic time is the corresponding transport time scale and the charac- 
teristic size of the wavetrains results from the saturation of the small nonlinear terms. The 
small amplitude slow envelope assumption, which is the key assumption that allow us to de- 
rive both systems of equations, forces the dispersive terms to be always small as compared 
with the advection terms (in other words, the nonzero group velocity turns this system into 
a transport dominated one) and therefore (with the scaling mentioned above) the NLCMEd 
must be considered only in the physically relevant regime e ^ 0. 

2. The NLCMEd are asymptotically nonuniform, in the sense that the resulting asymptotic 
model, in the e ^ limit, still contains the small parameter e. This is due to the fact 
that the NLCMEd include simultaneously two balances with different asymptotic order: 
one induced by the dominant transport terms and the other associated with the underlying 
effect of dispersion. This kind of asymptotically nonuniform amplitude equations have been 
previously derived in the context of water waves ||9|] and for the onset of the oscillatory 
instability in spatially extended dissipative systems 



HQ. 



3. Two spatial scales are present in the NLCMEd: transport scales (5xtrans ~ 1, and dispersive 
scales 5a:disp ~ \/W\ ^ 1 • The dispersive scales are small as compared with the transport 
scales but still large as compared with the wavelength of the basic resonant wavetrains in 
expression which, in the scaling we are using, is of the order of <^ 1, and therefore 
the slow envelope assumption is not violated. 
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4. If only transport scales are present, then the dispersion terms in the NLCMEd 
I^A^;!;!, I^-Bxa;! ~ ^ 1 can be safely neglected (they produce only a small quantita- 
tive correction that vanishes as e — 0) and the evolution of the system is well represented 
by the NLCME. On the other hand, if the small dispersive scales do develop, then the NL- 
CME do not correctly predict the dynamics system and the NLCMEd must be used instead. 
The onset of the dispersive scales is not a higher order, longer time effect; it takes place 
in the same timescale of the NLCME no matter how small the dispersion coefficient e is. 
Once the dispersive scales appear they typically spread all over the domain giving rise to 
very complicated spatio-temporal dynamics. This dispersive destabilization can be simply 
regarded as the standard modulational instability of the NLS-like dynamics that lays beneath 
the dominant transport induced dynamics. 

5. The stability of the family of uniform modulus solutions, known as continuous waves (CW), 
is drastically affected by dispersion. The stability predictions for the CW from the NLCME 
differ completely from those obtained from the NLCMEd for both signs of the dispersion 
coefficient e, no matter how small it may be. 

Despite of the results presented in (jsj] it appears that the NLCME continue to be used as the 
amplitude equations for the description of light propagation in FBG and for the weakly nonlinear 
dynamics of BEC in optical lattices without paying any attention to the effect of dispersion. In 
order to make clear that the correct amplitude equations are the NLCMEd, we have decided to 
carry out some numerical integrations of the full ID Maxwell- Lorenz equations (MLE) in a long 
fiber Bragg grating and check that the stability predictions for the CW given by the NLCME are 
wrong and that the dispersive NLCMEd give the correct results. 

This paper is organized as follows: in the following section we derive the explicit expressions 
of the coefficients of the NLCMEd from the MLE and, in the next and final section of this paper, 
we present some numerical integrations of the MLE starting from a perturbed CW and compare 
them with the CW stabiUty characteristics predicted by the NLCME and the NLCMEd. 



II. NLCMED DERIVATION FROM THE MLE 

Our formulation follows closely that of ref. [S]. We have decided to include here a quited 
detailed derivation of the NLCMEd from the MLE because we use rather new derivation procedure 
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that has the advantage of not requiring to assume any a priori relation among the different small 
parameters of the problem. 

We describe the propagation of light in a fiber with a periodic grating and a cubic nonlinearity 
using the one-dimensional Maxwell's equations 11121. 11311 for the evolution of the electromagnetic 



fields together with an anharmonic Lorentz oscillator model for the polarization (see e.g. [|5 
and references therein) 

dB _ dE 

dt dx ' 
dD _ dB 

dt dx ' 

D = eoE + P, 



P dt^ 



[1 - 2 An cos(27rx/Ag))P - -fP^ = eoxE. 



(6) 

(7) 
(8) 
(9) 



In the system above, the electric field E, the magnetic field B, the dielectric displacement D and 
the polarization P are scalar fields that depend on the spatial variable x and on time t. fiQ and 
eo denote, respectively, the permeability and the permittivity of the vacuum. The characteristic 
frequency fip accounts for the non instantaneous polarization response of the media. An and 
represent the strength and the period of the grating, that is, the strength and the period of the spatial 
periodic variation of the refractive index of the fiber (An measures the size of the nonuniformities 
of the refraction index relative to its mean value no, see Fig.[T]), x is the linear polarizability of the 
medium (n^ = 1 + x) 7 > is the coefficient of the nonlinear Kerr effect. 



Figure 1 : One dimensional fiber with a periodic variation of the refractive index. 

In order to simplify subsequent calculations it is convenient to make the system ©-(HI) nondi- 
mensional using the following rescalings: 



B = ^i2o/{eoi)B, D = {1/^)D, E={l/y^)E, P={1/^)P, 

X = {\g/Tx)x, t = {Xg/c7r)i, 



where (? = l/(eo/io) is the vacuum speed of light. After dropping tildes and eliminating D and 
B, the nondimensional MLE can be written in the form 

d^{E + P) d^E 



(10) 



- 2Ancos(2x))P + - \)E + ujIP\ (11) 

where the grating period is now equal to vr and the dimensionless finite time polarization response 
frequency is given by cj^ = f2pA^/(c^7r^). 

In the absence of grating, the linear propagation characteristics of a wavetrain of the form 



\E{x,t) 1 


hi 











are given by the following dispersion relation 





H 











u;t-ulie + ujlnl)+coie = 0, (13) 
which, for nl > 1, has four real roots of the form 

uJk = ±^{k' + ujy,)/2 ± y/(P + ooy.y/A - ^Ik^ (14) 
and associated eigenvectors 

(15) 

The four branches of the dispersion relation (fT4l) are plotted in Fig. [2l There are two different 
behaviors for large wavenumbers: one is dominated by the finite time polarization response of 
the medium, ujk — ^ ±^p as A; — ±oo, and the other, ujk — ^ as A; — > ±oo, corresponds to 
propagation like in the vacuum, without polarization effects. 

The small nonuniformities of the refractive index. An <^ 1, and the effect of the small nonlin- 
earity can be accounted for by allowing the wavetrains that resonate with the grating to be slowly 
modulated in space and time 

K)(A+(x,t)e"+^"* + A-(x,t)e-^^+*^*) + c.c. + . . . , (16) 
where 




Vo= { and a;=^/(l + ..>^)/2±V(l + ^X)74-^,^ (17) 



6 







____ . ^ 


, 



















k 



k 



k 



Figure 2: Sketch of the dispersion relation (fT4l) . 
and the weakly nonlinear level of this approach requires essentially that 
• ■ ■ < -C < < 1, 



< |Af I < \A^\ < 1 and An < 1, 



(18) 



that is, small amplitudes that depend slowly on space and time and small grating strength. The 
solution of eqs. (fTOl)-(fTTI) and the amplitude equations can be expanded in powers of the small 



quantities An, A^, A^, A^^,. . . as 




Vo{A+e 



ix+iujt 



+ A'e 



-ix+iuit 



) + C.C. + 



At = a^A+ + at At + at At, + ..., 



A: 



a^A +a^A,+a2A,, + ..., 



(19) 
(20) 
(21) 



which, once inserted into eqs. (fT0l)- (fTT1) . provide a linear nonhomogeneous system for the contri- 
bution of each order. For the resonant terms, i.e., those proportional to e='="='=*'^*, a condition must 
be satisfied to ensure that there are not secular terms in the short time scale. In other words, the 
linear problems corresponding to the resonant terms are singular and hence a solvability condition 
must be satisfied by the nonhomogeneous part; these solvability conditions yield the coefficients 
of the amplitude equations. 

Notice that only the resonant terms contribute to the amplitude equations and only the amplitude 
equation for A~^ has to be calculated because the corresponding equation for A" can be obtained 
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by simply applying the symmetry 



X 



-X 



A- 



(22) 



which comes from the spatial reflection symmetry of the original problem (fTOl)-(fTTI). 

The linear terms in the amplitude equations can be easily anticipated because they correspond 
to the Taylor expansion of the dispersion relation (fT4l) at A; = 1 (see e.g. lilSn ). 



duJk 



dk 



. 1 d'^LOk 



k=l 



k=l 



2 dk"^ 

The first coefficient obviously vanishes (see eq. (fTT]) ) and the second and third coefficients are, 
respectively, the group velocity and the higher order dispersion, which, after making use of eq. 
(fT3l) . can be written as 



duJk 



dk 



k=l 

'2, 



, A / t2 ■ 
(jj -Up 



id 



1 d^uJk 



-i- 



2 dk"^ 



Iuj\uj^-1){uj^-ujI){?,ujI + uo^) 

-i 

2 



(23) 
(24) 



which correspond to the group velocity and dispersion of the fiber without grating. 

The first order, resonant contributions of the grating to the expansion of the solution (fT6l) and 
to the amplitude equation (|20|) are of the form 



WAnA'e 



ix+iLdt 



and 



wAnA 



where the two component vector W is given by the following linear, singular nonhomogeneous 
problem 

f^2_i Lj^l ^Tnnl Til 

Vo. 

This system can be solved only if the right hand side is orthogonal to the solution of the adjoint 
problem 





W = -u;l 








Vo + 2iLJw 


1 


1 









1 







1 



ya 



and this solvability condition gives the value of the coefficient of the amplitude equation 

Ujil — UJ^) 2 



(25) 



The first order contributions of the nonlinear term. 



Ui A+\A+\ 2e^^+*^* , U2A+\A-\ 2e*^+^'^* 



and 



U2A+\A- 
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are computed similarly: the following linear problems are obtained for the vectors Ui and U2 










1 1 






+ 2iujui 










1 












1 1 




f/2 = -6^J 




+ 2iuJU2 










1 



V,. 



Vn. 



and, after applying the solvability condition, the resulting amplitude equation coefficients are given 
by 



Ui 



' —Up and U2 



.3tu(l-w2)3 2 



(26) 



The ratio U2 = 2uy could have been advanced; it is a well known result of the cubic nonlinearity 
of the problem [lisll . 

Collecting the coefficients above (|23l)-(|26l) and applying the spatial reflection symmetry (|22]) 
the resulting amplitude equations can be written as 



Af = VgA+ + idA+^ + wAnA- + + U2\A~\'^) + ... 



A: 



U2\A' 



+ .... 



(27) 
(28) 



It is important to emphasize that no particular scaling among the small size of the amplitudes, 
the small grating depth, the slow time and the large spatial scale has been used; only the slow 
envelope, weakly nonlinear assumption expressed in (fTSi) is actually required to obtain the above 
amplitude equations. 

We will consider the simplest possible geometrical configuration: propagation of light in a fiber 
ring with length L ^ 1. The spatial periodicity condition implies that the boundary conditions for 
y4+ and A~ are (see eq. (fT6l )) 



v4+(x + L)e^^ = A+{x,t), A-{x + L)e-''^ = A-{x,t). 



(29) 



Here 9 = L (mod27r) measures the mismatch between the natural wavelength of the resonant 
wavetrains (=2%) and the period of the domain, but we will confine ourselves to the particular case 
9 = 0, i.e., ring length equals to an integer multiple of the period of the wavetrains. 

There are two possible choices uj^ depending on the sign selected in (fTTI) . see Fig.[3l The group 
velocity (|23l) is positive in both cases (it is the slope of the curve cj^ at A; = 1 in Fig. |3]), but the 
sign of the dispersion coefficient d (|24|) . which is related to the curvature of the curve Uk in Fig. [3l 
changes. On the other hand, the nonlinear and grating terms have imaginary parts that are always 



negative; see eqs. (1251) and (l26l) and Fig. [3l and recall that, using the dispersion relation eq. (fT3l) 
for k = 1, the denominator can be written as lu"^ — cUp = uj'^[{u'^ ~ 1) + (^^ ~ '^p'^o)]- 

















Up ^ . ' 









k 



k = 1 



Figure 3: Detail of the dispersion relation ([Mi l with the two frequencies for = 1. 

In order to make the nonlinear and grating coefficients positive, we will work with the complex 
conjugates of the amplitudes and, to absorb some parameters of the problem, we will also perform 
the following rescalings 

x = Lx, t= {L/vg)i, = ^Jvg/{L\u2\)A^, (30) 

that, after dropping tildes, yield the scaled NLCMEd 

At= A+ +ieA+^ + iKA~ +iA+{a\A+\^ + \A~\^), (31) 
A- = -A- + ieA'^ + inA-^ + + (32) 

A^{x + l,t) = A^{x,t), (33) 

where e = —d/ (Lvg) ^ 1 is positive (negative) for u = uj~^(uj = uj~), the scaled grating strength 
K, = AnL\w\/vg ~ 1 is always positive, and the nonlinear coefficient cr = | (the standard NLCME 
are obtained by just by setting e = in the system above). 

III. NUMERICAL RESULTS 

In order to confirm that the correct stability predictions for the MLE are those given by the 
NLCMEd, we numerically integrate the complete MLE (fTOl)- (fTT]) in a large ring shaped fiber 
grating, that is, with periodic boundary conditions, 

E{x + L,t) = E{x, t), P{x + L,t) = P(x, t), 
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and L ^ 1. The MLE are integrated numerically as a system of four first order equations, using 



Fourier series in space and a 4th order Runge-Kutta scheme nl6\ for the time integration of the 



resulting ODEs. The linear diagonal terms are integrated implicitly and the nonlinear terms are 



computed in physical space using the 2/3 rule to remove the aliasing terms nllu . The number of 
modes used in the simulations presented is Afpourier = 1024 and the time step At = .01, and the 
Fourier transforms were performed using the FFTW routines jisl. 
The initial condition for all simulations is a CW 

= pCOSe 6'"*+'"^^ , 4" = p sin ^ giat+imx ^ 
K 1+0" K 1 — cr „ 

sm2t^ 2 sm2t/ 2 

where p > is the light intensity in the fiber and 9 g] — |, 0[U]0, | [ measures the ratio between the 
two counterpropagating wavetrains (see Is]), with a small superimposed perturbation. Once a CW 
has been selected (k, p and fixed) and the three MLE parameters cUp, and L are prescribed, 
the initial condition for the MLE is obtained from 











\-\ 





L\U2\ 

and the remaining MLE coefficient. An, and the dispersion coefficient of the NLCMEd, e, are 
given by 

An = K and e = (35) 

L\W\ LVg 

which can be computed after making use of (fTTI) , (|23l) , (|24l) . (|25l ) and (|26l ). 

We consider only two configurations because the MLE numerical integrations are rather CPU 
costly (large system length and very long final integration time). 

CASE 1 The initial CW parameters are k = 1, 6* = -f and p^ = 1. The NLCMEd results pre- 
sented in ref. [0] indicate that this CW is stable for negative dispersion and dispersively unstable 
for positive dispersion, while, according to the NLCME, this CW is always stable. The numer- 
ical integrations of the MLE presented in Fig. |4] correspond to the parameters uj"^ = 1, = 2, 
L = 128n (i.e., there are 128 grating oscillations inside the fiber ring). The first and second plot 
correspond, respectively, to cu = u~ and lu = (that is, to negative and positive e in the NL- 
CMEd (see eq (1351) and Fig.. [3])) with the MLE grating strength. An, that results from eq. (I35l) . 
They show the time evolution of the spatial norm of the electric field. 
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and look like a solid black patch due to the fact that \ \E\ \ oscillates very fast in time. In agreement 
with the NLCMEd predictions, the CW is stable for negative dispersion (first plot in Fig. and 
unstable for positive dispersion (the instability growth can be appreciated from t = 60000 on in 
the second plot of Fig. U). The corresponding spatial profiles ofEatt = 75000 are given in the 
third and fourth plots of Fig.Hl for negative dispersion (third plot) a perfectly uniform amplitude 
oscillatory pattern is obtained (the CW pattern) but, for positive dispersion, a modulation is clearly 
present (fourth plot). In order to be sure that this is a dispersive instability we have repeated the 
unstable MLE simulation in a four times longer domain (L = 256(27r)). The resulting spatial 
profile of -E at t = 160000 is shown in the last plot of Fig. IH Notice how the number of basic 
wavelengths is now four times higher but the number of wavelengths of the modulation only ap- 
proximately doubles (increases from 5 to 9), confirming the dispersive character of the instability 



whose characteristic size scale as v -^^ (see ref. [|8l]). 

CASE 2 The CW parameters are now k = 1, 9 = ^ and = 1, and the MLE parameters 
are the same as in the above case: cUp = 1, = 2 and L = 1287r. The first and third plot of 
Fig. [5] correspond to positive dispersion and indicate that the CW is now stable. The dispersion 
is negative in the second and fourth plot where the destabilization of the CW can be clearly seen 
both in the time evolution of ll-Ell and in the dispersive modulations that the spatial profile of 
E displays. This is again in perfect agreement with the linear stability results obtained from the 
NLCMEd iSQ (the NLCME again wrongly labeled this CW as always stable). 

In conclusion, the numerical simulations of the MLE indicate that the NLCME fail to describe 
the system evolution if dispersive instabilities (that cannot be detected using the NLCME formu- 
lation) come into play. In this case the higher order dispersion effects must be taken into account, 
and the amplitude equations that do correctly predict the weakly nonlinear dynamics of light prop- 
agation in FBG are the asymptotically nonuniform NLCMEd dH)-®. 
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Figure 4: MLE simulation results starting from a CW (k = 1, 9 = —j and = 1) with a 10 ^ perturbation. 
From top to bottom: time evolution of the spatial norm of E for lu~ and u>~^, spatial profiles of E at 
t = 75000 for uj~ and lo+, and spatial profile of ^ at t = 160000 for and L = bUir. 
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Figure 5: MLE simulation results starting from a CW (k = 1, 9 = j and = 1) with a 10~^ perturbation. 
From top to bottom: time evolution of the spatial norm of E for and lo~ , spatial profile of £^ at t = 75000 
and io~^, and spatial profile of -E at t = 30000 for a;~. 
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